In [1]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
In [2]:
from gevi_classes import *

/Users/GP1514/.pyenv/versions/anaconda3-2.5.0/lib/python3.5/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)
In [3]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 1000;
In [4]:
# instantiate utility class
gr = Graph()

Load DATA

Data collected in 5 X 1min, with 1min pause in between 2921 images obtained in 60s

In [5]:
username = os.path.expanduser('~').split('/')[-1]
if username == "GP1514":
    print("At Imperial")
    mouseAPath = '/Volumes/DATA/GEVI/MouseA/'
    mouseBPath = '/Volumes/DATA/GEVI/MouseB/'
    mouseM1217FPath = '/Volumes/DATA/GEVI/M1217F/'
    mouseM1223MPath = '/Volumes/DATA/GEVI/M1223M/'
else:
    print("Using laptop")
    mouseAPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 27/'
    mouseBPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 22/'
    mouseCPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 28/'

At Imperial
In [6]:
discard = loadDic()
In [7]:
discard
Out[7]:
{'MouseA': {3: [0, 1, 3], 6: [1]},
 'MouseB': {3: [1, 4]},
 'MouseC': {},
 'MouseM1217F': {},
 'MouseM1223M': {}}
In [8]:
mouseA = Mouse('mouseA', mouseAPath, [3,4,5,6],discard['MouseA'] )
# mouseA = Mouse('mouseA', mouseAPath, [3],discard['MouseA'] )
mouseB = Mouse('mouseB', mouseBPath, [2,3,4,5],discard['MouseB'] )

mouseM1217F = Mouse('mouseM1217F', mouseM1217FPath, range(1,14),discard['MouseC'] )
mouseM1223M = Mouse('mouseM1223M', mouseM1223MPath, range(21,36),discard['MouseC'] )
In [9]:
# mouseA.loadData()
# mouseB.loadData()
mouseM1217F.loadData()
mouseM1223M.loadData()
Out[9]:
1
In [10]:
# os.chdir(mouseM1217F.path)
# dataMask = Tree()
# for file in glob.glob("*traceMask.hd5f"):
#     print(file)
#     dataMask[int(file[3:6])][int(file[12:15])] = tables.openFile(mouseM1217F.path + file, mode='r')
In [11]:
# dataMask[1][1].root.mask
mouseM1217F.mouseMask[0,100]
Out[11]:
1
In [12]:
mask = 1-mouseM1217F.mouseMask[:]
x = mouseM1217F.experiments[0].repeats[0].ratio
xm = [ma.masked_array(x[i],np.swapaxes(mask,0,1)) for i in range(x.shape[0])]
In [13]:
# xm[0].shape
res = list(map(np.nanmean,xm))
# plt.plot(res)
In [14]:
plt.plot(res)
rm = mouseM1217F.experiments[0].repeats[0].mRatio
plt.plot(rm)
plt.xlim([0,100])
Out[14]:
(0, 100)
In [15]:
# mouseB.experiments[1].repeats[0].setInfo('End of transition, then discard')
# mouseB.experiments[2].repeats[0].setInfo('Transition to anesthesia')
# mouseA.experiments[0].repeats[2].setInfo('Keep')
# mouseA.experiments[2].repeats[1].setInfo('Transition to desynchronization: short periods of silence trigger recovery of hemodynamics')
# mouseA.experiments[3].repeats[1].setInfo('Woke up - Discard')
# mouseA.experiments[3].repeats[2].setInfo('Discard')
# mouseA.experiments[3].repeats[3].setInfo('Discard')
# mouseA.experiments[3].repeats[4].setInfo('Discard')
# mouseA.getInfo()
# mouseB.getInfo()
In [16]:
#len(mouseM1217F.experiments[0].repeats[0].mRatio)
len(mouseM1217F.experiments[0].repeats[0].ratio)
Out[16]:
11800

Plot Data

Mouse M1217F

In [17]:
gr.plotHV(mouseM1217F,mask=True)

Mouse M1223M

In [18]:
gr.plotHV(mouseM1223M, mask=True)

Transfer functions

$V * \alpha = H$

$\mathscr{F}$ transform : $\hat{V}.\hat{\alpha} = \hat{H}$

$\alpha = \mathscr{F}^{-1} \frac{\hat{H}}{\hat{V}}$

ALPHA CURVE FITTING

$\alpha(t)=\frac{\tau_0}{t+0.1} - \frac{t}{\tau_1}*e^{\tau_3-\frac{t}{\tau_2}}$

Extra functions

In [19]:
def meanCorH(mouse):
    corH = np.mean(flatten(mouse.experiments.corH))
    return corH

def meanCorR(mouse):
    corR = np.mean(flatten(mouse.experiments.corR))
    return corR
    
# optimization for mouseA
def fn(p_est):
    alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
    gr.computeCorr(mouseM1223M, alpha)
    corR = np.mean(flatten(mouseM1223M.experiments.corR))
    res = (1 - corR)**2
#     print('%.2f \t\t t0:%.5f \t t1:%.5f \t t2:%.5f \t t3:%.5f'%(res, p_est[0],p_est[1],p_est[2], p_est[3]))
    return res

def fnH(p_est):
    alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
    gr.computeCorr(mouseM1223M, alpha)
    corH = np.mean(flatten(mouseM1223M.experiments.corH))
    res = (1 - corH)**2
#     print('%.2f \t\t t0:%.5f \t t1:%.5f \t t2:%.5f \t t3:%.5f'%(res, p_est[0],p_est[1],p_est[2], p_est[3]))
    return res

def computeAndPlotCorr(mouse, alpha):
        for i,exp in enumerate(mouse.experiments):
            for j,rep in enumerate(exp.repeats):
                gr.plotRModel(mouse, exp = i, rep = j,alpha = alpha )
                gr.plotHModel(mouse, exp = i, rep = j,alpha = alpha )
                gr.corrCoeff(mouse, exp = i, rep = j,alpha = alpha )
                
x = np.real(xax(mouseM1217F.experiments[0].repeats[0].mRatio, 60000*4))
In [20]:
# r= mouseA.experiments[0].repeats[0].ratio
# out=np.ones((60,60))
In [21]:
# filtering parameters
mouseM1217F.minFreqAlpha = 1
mouseM1217F.maxFreqAlpha = 4000
mouseM1223M.minFreqAlpha = 1
mouseM1223M.maxFreqAlpha = 4000
# mouseA.window = 10
gr.fHmax = 120
gr.fRmax = 800

Tranfer functions mouse M1217F

In [22]:
gr.plotTF(mouseM1217F, mask=True)
Out[22]:
1
In [23]:
gr.plotmTF(mouseM1217F, mask=True)
Out[23]:
1

Transfer functions mouse M1223M

In [24]:
gr.plotTF(mouseM1223M, mask=True)
Out[24]:
1
In [25]:
gr.plotmTF(mouseM1223M, mask=True)
Out[25]:
1

Voltage and Hemo reconstruction : with mean alpha function

/!\ Models from the same mouse, different experiments: Model from M1217F, exp11, applied on M1217F, exp3

In [26]:
# alpha = mouseA.experiments[3].repeats[3].meanAlphaModel
# mouseM1217F.experiments[0].repeats[0].getAlpha()
alpha = mouseM1217F.experiments[10].repeats[0].alphaMasked
gr.plotRModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )

(839.678011501-126.396537585j)

Model from M1223M, exp1, applied on M1217F, exp3

In [27]:
alpha = mouseM1223M.experiments[0].repeats[0].meanAlphaModel
gr.plotRModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )

(1808.03680337-310.357729984j)

Voltage and Hemo reconstruction : with mean alpha model function

Model from M1223M, exp5, applied on M1217F, exp3

In [28]:
# alpha = mouseA.experiments[3].repeats[3].meanAlphaModel
# gr.fHmax=50
# gr.fRmax=50
# mouseM1223M.experiments[2].repeats[0].getAlpha()
alpha = mouseM1223M.experiments[4].repeats[0].meanAlphaModel
gr.plotRModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha, window=10 )
gr.plotHModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha, window=10 )
plt.figure()
plt.plot(alpha)

(2648.5033912-375.71570331j)
Out[28]:
[<matplotlib.lines.Line2D at 0x2432b65f8>]

Distribution of alpha model parameters $\tau_0, \tau_1, \tau_2, \tau_3$

$v * \alpha = h$

$\alpha(t)=\frac{\tau_0}{t+0.1} - \frac{t}{\tau_1}*e^{-\frac{t}{\tau_2}}$

In [29]:
# alpha = mouseB.experiments[2].meanAlphaModel
# fontsize=16
# plt.plot(xax(alpha,60000),alpha)
# plt.xlabel('Time [s]')
# plt.title('Alpha model of mean TF mouse B')
# plt.text(30,-0.01, r'$ \tau_{0} = %.2g$' %mouseB.experiments[2].meanAlphaParams[0], fontsize =fontsize)
# plt.text(30,-0.02, r'$ \tau_{1} = %.2g$' %mouseB.experiments[2].meanAlphaParams[1], fontsize =fontsize)
# plt.text(30,-0.03, r'$ \tau_{2} = %.2g$' %mouseB.experiments[2].meanAlphaParams[2], fontsize =fontsize)
In [30]:
# gr.computeCorr(mouseM1217F, alpha=None)
# gr.computeCorr(mouseM1223M, alpha=None)
# gr.computeAlphaModels(mouseM1217F)
# gr.computeAlphaModels(mouseM1223M)
# gr.plotParamsIndex(mouseM1217F)
# gr.plotParamsIndex(mouseM1223M)
# gr.plotParamsCor(mouseM1217F)
# gr.plotParamsCor(mouseM1223M)
# print([meanCorR(mouseM1217F),meanCorH(mouseM1217F),meanCorR(mouseM1223M), meanCorH(mouseM1223M)])

Optimize for hemo

In [31]:
resultHemo = minimize(fnH, x0=[1,1,1], method = 'COBYLA',
                      options={'gtol': 1e-8, 'disp': True})
print(resultHemo)

     fun: 0.35133626601258383
   maxcv: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 143
  status: 1
 success: True
       x: array([ 0.67578356,  0.58751443,  0.98884867])
In [32]:
p_est = resultHemo.x
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
plt.plot(xax(alpha,4*60000),alpha)
# alpha2 = mouseB.experiments[2].meanAlphaModel
# plt.plot(xax(alpha2,60000),alpha2)
Out[32]:
[<matplotlib.lines.Line2D at 0x204fcda58>]
In [33]:
# computeAndPlotCorr(mouseA, alpha)
gr.plotRModel(mouseM1217F,2,0,alpha)
gr.plotHModel(mouseM1217F,2,0,alpha)
# gr.plotHModel(mouseA,2,4,mouseB.experiments[2].meanAlphaModel)
# gr.plotRModel(mouseA,2,4,mouseB.experiments[2].meanAlphaModel)
def printl(tab):
    for t in tab:
        print(t)

(2163.07704027-643.836581104j)

Optimize for voltage

In [34]:
# bounds = [(-10,10),(-10,10),(-10,10),(-10,10)]
# resultVolt = differential_evolution(fn, bounds, init = 'random')
resultVolt = minimize(fn, method = 'COBYLA', x0=[-1e-2,1,10,1],
                      options={'gtol': 1e-6, 'disp': True})
print(resultVolt)

     fun: 0.66311394545387858
   maxcv: 0.0
 message: 'Maximum number of function evaluations has been exceeded.'
    nfev: 1000
  status: 2
 success: False
       x: array([  0.11048601,   0.61678145,  13.58790947,   1.04188789])
In [35]:
p_est = resultVolt.x
print(p_est)
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]

[  0.11048601   0.61678145  13.58790947   1.04188789]
In [36]:
fontsize=10
plt.title('Optimized Alpha model of mean TF mouse B, exp %d'%(2+mouseB.start))
plt.text(30,0.5, r'$ \tau_{0} = %.2g$' %p_est[0], fontsize =fontsize)
plt.text(30,1.5, r'$ \tau_{1} = %.2g$' %p_est[1], fontsize =fontsize)
plt.text(30,2.5, r'$ \tau_{2} = %.2g$' %p_est[2], fontsize =fontsize)
plt.plot(xax(alpha,4*60000),alpha)
Out[36]:
[<matplotlib.lines.Line2D at 0x247d31f28>]
In [37]:
gr.plotRModel(mouseM1217F,2,0,alpha)
gr.plotHModel(mouseM1217F,2,0,alpha)
##
gr.plotRModel(mouseM1217F,9,0,alpha)
gr.plotHModel(mouseM1217F,9,0,alpha)

(317.107582274-874.801718299j)
(1242.4256318-139.394226031j)

With the last model developed for mice A and B (from optimizing voltage model on voltage data from mouse B)

In [38]:
p_est = [0.4121683, 1.09384498, 3.68435018,  0.98029149]
In [39]:
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
gr.plotRModel(mouseM1217F,2,0,alpha)
gr.plotHModel(mouseM1217F,2,0,alpha)

(-84.5694100206-1302.90938144j)
In [40]:
gr.plotRModel(mouseM1217F,9,0,alpha)
gr.plotHModel(mouseM1217F,9,0,alpha)

(1106.72936772-324.781800756j)
In [41]:
%timeit fn([1,2,3])

1 loop, best of 3: 254 ms per loop

TF on different experiments: M1217, exp 10

Model from M1223M, exp1, applied on M1217F, exp10

In [42]:
alpha = mouseM1223M.experiments[0].repeats[0].meanAlphaModel
plt.figure(figsize=(10,1))
plt.plot(alpha)
gr.plotRModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )

(1502.53001311-1028.47669135j)
In [43]:
alpha = mouseM1223M.experiments[4].repeats[0].meanAlphaModel
plt.figure(figsize=(10,1))
plt.plot(alpha)
gr.plotRModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )

(1905.59979673-85.6129856895j)
In [ ]: